/// 网格信息文件
#ifndef __cellinf
#define __cellinf

#include <iostream>
#include <vector>
#include <string>
#include "type.h"
#include "particles.h"
#include <cstring>
#include "physicalconstants.h"
#include "stdio.h"
#include "events.h"
#include "common.h"
#include "userdef.h"
#include "denfunc.h"
#ifdef __OMP
#include "omp.h"
#endif
using namespace std;
///用来存储MPI传输的Ghost cell数据，只有场，没有粒子
struct __Cdata{
    /// 0 左下， 1, 左上， 2 右下， 3, 右上
    //double data_rho[4];
    /// 0 先jz，后jx[0], jx[1], jy[0], jy[1]
    //double data_j[5];
    /// 电磁场同例
    double data_e[5];
    double data_b[5];
};

double get_density(const __Vect2<double>& position, const int& sp);
///单个cell的场信息
class __Cellfield{
public:
    /// 下标规则，总是先下后上，先左后右
    __Vect2<double> exfield;
    __Vect2<double> bxfield;
    //__Vect2<double> jxfield;
    __Vect2<double> eyfield;
    __Vect2<double> byfield;
    //__Vect2<double> jyfield;
    /// 左和下
    //__Vect2<double> rhofieldl;
    /// 右和上
    //__Vect2<double> rhofieldr;
    double ezfield;
    double bzfield;
    //double jzfield;
    __Cellfield(){
        __Vect2<double> zero(0, 0);
        exfield = zero;
        eyfield = zero;
        bxfield = zero;
        byfield = zero;
        //jxfield = zero;
        //jyfield = zero;
        ezfield = 0;
        bzfield = 0;
        //jzfield = 0;
        //rhofieldl = zero;
        //rhofieldr = zero;
    }
    __Cellfield(const __Cellfield& other){
        exfield = other.exfield;
        eyfield = other.eyfield;
        ezfield = other.ezfield;
        bxfield = other.bxfield;
        byfield = other.byfield;
        bzfield = other.bzfield;
        /***
        jxfield = other.jxfield;
        jyfield = other.jyfield;
        jzfield = other.jzfield;
        rhofieldl = other.rhofieldl;
        rhofieldr = other.rhofieldr;
        ***/
    }
    __Cellfield(__Cdata& cdata);
    ~__Cellfield(){}
    __Cellfield operator=(const __Cellfield& other){
        this->exfield = other.exfield;
        this->eyfield = other.eyfield;
        this->ezfield = other.ezfield;
        this->bxfield = other.bxfield;
        this->byfield = other.byfield;
        this->bzfield = other.bzfield;
        /***
        this->jxfield = other.jxfield;
        this->jyfield = other.jyfield;
        this->jzfield = other.jzfield;
        this->rhofieldl = other.rhofieldl;
        this->rhofieldr = other.rhofieldr;
        ***/
        return *this;
    }
    ///所有数据清空为 0
    void clear();
};

void __Cellfield::clear(){
    __Vect2<double> zero(0, 0);
    exfield = zero;
    eyfield = zero;
    bxfield = zero;
    byfield = zero;
    //jxfield = zero;
    //jyfield = zero;
    ezfield = 0;
    bzfield = 0;
    /***
    jzfield = 0;
    rhofieldr = zero;
    rhofieldl = zero;
    ***/
    }

__Cellfield::__Cellfield(__Cdata& cdata)
{
    /***
    rhofieldl.member[0] = cdata.data_rho[0];
    rhofieldl.member[1] = cdata.data_rho[1];
    rhofieldr.member[0] = cdata.data_rho[2];
    rhofieldr.member[1] = cdata.data_rho[3];

    jzfield = cdata.data_j[0];
    ***/
    ezfield = cdata.data_e[0];
    bzfield = cdata.data_b[0];

    /**
    jxfield.member[0] = cdata.data_j[1];
    jxfield.member[1] = cdata.data_j[2];
    jyfield.member[0] = cdata.data_j[3];
    jyfield.member[1] = cdata.data_j[4];
    ***/

    exfield.member[0] = cdata.data_e[1];
    exfield.member[1] = cdata.data_e[2];
    eyfield.member[0] = cdata.data_e[3];
    eyfield.member[1] = cdata.data_e[4];


    bxfield.member[0] = cdata.data_b[1];
    bxfield.member[1] = cdata.data_b[2];
    byfield.member[0] = cdata.data_b[3];
    byfield.member[1] = cdata.data_b[4];
}

/// 单个cell内部信息，包括场和粒子两部分
class __Cellinfo{
    public:
       __Cellinfo(){}
       /// 网格位置
       __Vect2<int> position;
       /// 网格左下角坐标
       __Vect2<double> xymin;
       /// 网格右下角坐标
       __Vect2<double> xymax;
       /// 网格中的粒子束
       vector<__Ptclslist> plist;
       /// 网格中的场
       __Cellfield fields;
       /// particle species total
       int species;
       /// 粒子种类启用标志
       /// 粒子种类密度
       vector<double>  density = {};
       /// area of this cell
       double area;
       /// 0:synchrotron, 1:bremsstrahlung, 2:bw_pair, 3:bh_pair
       int events_types[5];
       __Cellinfo(__Vect2<int>& position, __Vect2<double>& xymin, __Vect2<double>& xymax, \
               vector<double>& dens, const int& species):position(position), xymin(xymin),\
                                                           xymax(xymax), density(dens), species(species){}
       __Cellinfo(__Vect2<int>& position, __Vect2<double>& xymin, __Vect2<double>& xymax, \
               const int& species):position(position), xymin(xymin), xymax(xymax), species(species){
                                                                for(int i = 0; i < species; i ++)
                                                                {
                                                                    density.push_back(0);
                                                                }
       }
       /// 使用其它网格数据进行初始化
        __Cellinfo(const __Cellinfo& other){
           position = other.position;
           xymin = other.xymin;
           xymax = other.xymax;
           plist = other.plist;
           fields = other.fields;
           density = other.density;
           species = other.species;
       }
        /// 可以使用=来直接对网格赋值
       __Cellinfo operator=(const __Cellinfo& other){
           this->position = other.position;
           this->xymin = other.xymin;
           this->xymax = other.xymax;
           this->plist = other.plist;
           this->fields = other.fields;
           this->density = other.density;
           this->species = other.species;
           return *this;
       }

       /// 网格粒子推进
       void advance(const double& delta_t, vector<vector<__Ptcls> >& templist, const int& myrank);
       /// 网格sp种类例子，偏移offset处粒子，删除
       void pop(const int& sp, const int& offset);
       /// 将myptcl加入到当前网格中
       void push(const __Ptcls& myptcl);
       /// 压缩空间
       void squeeze();
       /// 将rho和j weight到grid上
       void weight_to_grid();
       /// weight 场到坐标，粒子推动使用
       __Vect3<double> weight_field_to_coords(const int& type, const __Vect2<double>& coords);
        // type 0--efield, type 1--bfield
        /// 设置当前网格的场值
       void set_field(const __Cellfield& otherfield);
       /// 初始化当前网格，假定density已经给定
       void init(const vector<int>& ppc);
       /// 初始化sp类粒子，给定ppc和dens
       void init(const vector<int>& ppc, const vector<Particles>& spptcl, const int& disable_em);
       /// 设这sp类粒子的速度
       void set_velocity(const int& sp, const __Vect3<double>& velocity);
       /// 设置当前网格的粒子种类密度，dens是数组
       void set_density(const vector<double>& dens);
       /// 将场插值到粒子，推动粒子使用
       void weight_field_to_ptcls();
       /// 同步辐射，采用光学深度算法
       void synchrotron(const double& delta_t, vector<vector<double> >& tables1d, \
               vector<vector<vector<double> > >& tables2d);
       /// 韧致辐射，采用MC
       void bremsstrahlung(const double& delta_t, \
               vector<vector<double> >& tables1d, \
               vector<vector<vector<double> > >& tables2d);
       /// 光子光子产生正负电子对Breit_wheeler过程，采用光深，以后可以修改为与背景温度相关的MC过程。
       void bw_pair(const double& delta_t, \
               vector<vector<double> >& tables1d, \
               vector<vector<vector<double> > >& tables2d);
       /// 光子在原子核的库仑场中产生正负电子对，采用MC
       void bh_pair(const double& delta_t, \
               vector<vector<double> >& tables1d, \
               vector<vector<vector<double> > >& tables2d);
       /// 粒子信息输出重载算符，可以直接输出当前粒子的信息
       friend ostream& operator<<(ostream& out, const __Cellinfo& other){
           out<<"position "<<other.position<<"\t"<<other.xymin<<"\t"<<other.xymax<<endl;
           // 注意这个地方的定义必须是引用
           for(int i = 0; i < other.species; i++){
                for(int j = 0; j < other.plist[i].member.size(); j++){
                    out<<other.plist[i].member[j]<<endl;
               }
           }
           out<<"________________________"<<endl;
           return out;
       }

       /// 计算各种粒子的密度(本网格内部)

       ~__Cellinfo(){}
};


void __Cellinfo::synchrotron(const double& delta_t, vector<vector<double> >& tables1d, \
               vector<vector<vector<double> > >& tables2d)
{
    //cout<<"size "<<plist[0].size()<<endl;
    /// 电子
    if(plist[0].member.size() > 0)
    {
        for(int i = 0; i < plist[0].member.size(); i ++)
        {
            __Ptcls child0;
            int status = synchrotron_rad(delta_t, plist[0].member[i], \
                    child0, tables1d, tables2d);
#ifndef PHOTON_DYNAMICS
        status = 0;
#endif
            if(status == 1)
            {
                plist[3].push(child0);
                //cout<<"new photon generated"<<endl;
            }
            else {
                child0.total_id -= 1;
            }
        }
    }
    if(plist[2].member.size() > 0)
    {
        //cout<<" starting positron "<<endl;
        /// 正电子
        for(int i = 0; i < plist[2].member.size(); i ++)
        {
            __Ptcls child0;
            int status = synchrotron_rad(delta_t, plist[2].member[i], \
                    child0, tables1d, tables2d);
#ifndef PHOTON_DYNAMICS
        status = 0;
#endif
            if(status == 1)
            {
                plist[3].push(child0);
                //cout<<" new photon generated by posiron"<<endl;
            }
            else {
                child0.total_id -= 1;
            }
        }
    }
}
void __Cellinfo::bremsstrahlung(const double& delta_t, \
        vector<vector<double> >& tables1d, \
        vector<vector<vector<double> > >& tables2d)
{
    for(int i = 0; i < plist[0].member.size(); i ++)
    {
        __Ptcls child;
        int status = bremsstrahlung_rad(delta_t, density[1], plist[0].member[i], \
                child, tables1d, tables2d);
#ifndef PHOTON_DYNAMICS
        status = 0;
#endif
        if(status == 1)
        {
            plist[3].push(child);
        }
        else {
            child.total_id -= 1;
        }
    }
    for(int i = 0; i < plist[2].member.size(); i ++)
    {
        __Ptcls child;
        int status = bremsstrahlung_rad(delta_t, density[1], plist[2].member[i], \
                child, tables1d, tables2d);
#ifndef PHOTON_DYNAMICS
        status = 0;
#endif
        if(status == 1)
        {
            plist[3].push(child);
        }
        else {
            child.total_id -= 1;
        }
    }
}
void __Cellinfo::bw_pair(const double& delta_t, vector<vector<double> >& tables1d, \
               vector<vector<vector<double> > >& tables2d)
{
    //cout<<" in bw pair "<<endl;
    ///　因为会删除粒子因此，必须倒序才能避免跳过粒子，或者对空list操作
    if(plist[3].member.size() > 0)
    {
        for(int i = plist[3].member.size() - 1; i >= 0; i --)
        {
            __Ptcls child0, child1;
            int status = bw_pair_create(delta_t, plist[3].member[i], \
                    child0, child1, tables1d, tables2d);
            if(status == 2)
            {
                ///产生粒子后这个list 会减小。
                plist[0].push(child0);
                plist[2].push(child1);
                plist[3].pop(i);
            }
            else {
                child0.total_id -= 2;
            }
        }
    }
}

void __Cellinfo::bh_pair(const double& delta_t, \
        vector<vector<double> >& tables1d, \
        vector<vector<vector<double> > >& tables2d)
{
    if(plist[3].member.size() > 0)
    {
        for(int i = plist[3].member.size() - 1; i >= 0; i --)
        {
            __Ptcls child0, child1;
            int status = bh_pair_create(delta_t, density[1], plist[3].member[i], \
                    child0, child1, tables1d, tables2d);
            if(status == 2)
            {
                plist[0].push(child0);
                plist[2].push(child1);
                plist[3].pop(i);
            }
            else {
                child0.total_id -= 2;
            }
        }
    }
}




void __Cellinfo::set_density(const vector<double>& dens)
{
    density = dens;
}


/// 一阶插值函数，这里采用面积插值, mycoords处的场插值到other，或者相反（结果是一样的），给定cell_length
double weightfunc(const __Vect2<double>& mycoords, const __Vect2<double>& other, \
        const __Vect2<double>& cell_length){
    __Vect2<double> temp(0, 0);
    temp = (mycoords - other).Abs();
    //cout<<"temp "<<temp<<endl;
       if(cell_length > temp){
           return (cell_length.member[0] - temp.member[0]) * \
               (cell_length.member[1] - temp.member[1])\
               / (cell_length.member[0] * cell_length.member[1]);
       }
    else {
        return 0;
    }
}


/// 二阶插值,没有使用，需要矫正

double weightfunc1d(const __Vect2<double>& pcoords, const __Vect2<double> ccoords, \
        const __Vect2<double>& cell_length)
{
    __Vect2<double> temp = pcoords - ccoords;
    double x, y;
    if(abs(temp.member[0]) < 0.5 * cell_length.member[0])
    {
        x = 0.75 - temp.member[0] * temp.member[0];
    }
    else if(abs(temp.member[0]) < 1.0 * cell_length.member[0])
    {
        x =  0.5 * (0.5 + abs(temp.member[0])) * (0.5 + abs(temp.member[0]));
    }
    else{
        return 0;
    }
    if(abs(temp.member[1]) < 0.5 * cell_length.member[1])
    {
        y = 0.75 - temp.member[1] * temp.member[1];
    }
    else if(abs(temp.member[1]) < 1.0 * cell_length.member[1])
    {
        y =  0.5 * (0.5 + abs(temp.member[1])) * (0.5 + abs(temp.member[1]));
    }
    else{
        return 0;
    }
    return x * y;
}





__Vect3<double> __Cellinfo::weight_field_to_coords(const int& type, const __Vect2<double>& coords){
    //cout<<" start type "<<type<<endl;
    __Vect3<double> retfields(0, 0, 0);
    __Vect2<double> dxy = xymax - xymin;
    __Vect2<double> ucenter(xymin.member[0] + dxy.member[0] / 2.0, xymax.member[1]);
    __Vect2<double> dcenter(xymin.member[0] + dxy.member[0] / 2.0, xymin.member[1]);
    __Vect2<double> lcenter(xymin.member[0], xymin.member[1] + dxy.member[1] / 2.0);
    __Vect2<double> rcenter(xymax.member[0], xymin.member[1] + dxy.member[1] / 2.0);
    __Vect2<double> center = xymin + dxy / 2.0;
    double wuc, wdc, wlc, wrc, wz;
    //cout<<" start weight func "<<endl;
    //cout<<" ucenter " <<ucenter<<endl<<" coords "<<coords<<endl<<" dxy "<<dxy<<endl;
    wuc = weightfunc(ucenter, coords, dxy);
    //cout<<" wuc over "<<endl;
    wdc = weightfunc(dcenter, coords, dxy);
    wlc = weightfunc(lcenter, coords, dxy);
    wrc = weightfunc(rcenter, coords, dxy);
    wz = weightfunc(center, coords, dxy);
    //cout<<" weight func finish "<<endl;
    switch(type)
    {
        // x: u/d center up, down
        // y: l/r center left, right
        // z: center

        case 0:
            retfields.member[0] += fields.exfield.member[0] * wdc + fields.exfield.member[1] * \
                                   wuc;
            retfields.member[1] += fields.eyfield.member[0] * wlc + fields.eyfield.member[1] * \
                                   wrc;
            retfields.member[2] += fields.ezfield * wz;
            //cout<<" e finish "<<endl;
            break;
        case 1:
            retfields.member[0] += fields.bxfield.member[0] * wdc + fields.bxfield.member[1] * \
                                   wuc;
            retfields.member[1] += fields.byfield.member[0] * wlc + fields.byfield.member[1] * \
                                   wrc;
            retfields.member[2] += fields.bzfield * wz;
            //cout<<" b finish "<<endl;
            break;
        default:
            break;
    }
    //cout<<"retfields is "<<retfields<<endl;

    return retfields;
}



void __Cellinfo::weight_field_to_ptcls(){
    for(int i = 0; i < species; i ++){
        for(int j = 0; j < plist[i].size(); j ++){
            //cout<<"list "<<i<<" of p "<<j<<" is weighting "<<endl;
            plist[i].member[j].efield = weight_field_to_coords(0, plist[i].member[j].position);
            //cout<<" coords e finishi "<<endl;
            plist[i].member[j].bfield = weight_field_to_coords(1, plist[i].member[j].position);
            //cout<<" coords to b finish "<<endl;
        }
    }
}


void __Cellinfo::init(const vector<int>& ppc, const vector<Particles>& spptcl, const int& disable_em){
    //cout<<"with xymin "<<xymin<<endl;
    __Vect2<double> zero2(0, 0);
    __Vect3<double> zero3(0, 0, 0);
    __Vect2<double> xydelta = xymax - xymin;
    __Vect2<double> xycenter = xymin + xydelta * 0.5;
    area = xydelta.member[1] * xydelta.member[0];
    //cout<<"xydelta "<<xydelta<<endl;
    __Cellfield zerofield;
    fields = zerofield;
    for(int i = 0; i < species; i++){
        //cout<<"init species "<<spptcl[i].name<<" with density "<<density.size()<<"\t"<<ppc.size()<<endl;
        __Ptclslist mylist(xydelta, i);
        if(ppc[i] > 0 && abs(density[i])>1e-30 && !disable_em){
            __Vect2<double> pos(0, 0);
            long weight = long(density[i] * area / ppc[i]);
            for(int j = 0; j < ppc[i]; j ++)
            {
                __Ptcls myptcl(spptcl[i]); 
                myptcl.weight = weight;
                //cout<<"pid "<<myptcl.total_id<<"\t sp "<<i<<endl;
                pos = xymin + xydelta * random_generator();
                myptcl.position = pos;
                myptcl.old_pos = pos;
                mylist.push(myptcl);
            }
        }
        plist.push_back(mylist);
    }
}

void __Cellinfo::set_velocity(const int& sp, const __Vect3<double>& velocity){
    for(int i = 0; i < plist[sp].size(); i ++){
        plist[sp].member[i].velocity = velocity;
    }
}



void __Cellinfo::squeeze(){
    for(int i = 0; i < plist.size(); i ++)
    {
        plist[i].squeeze();
    }
}


void __Cellinfo::push(const __Ptcls& myptcl)
{
    //cout<<"pushing into cell "<<position<<endl;
    for(int i = 0; i < plist.size(); i++)
    {
        //cout<<"trying sp "<<i<<endl;
        //cout<<"sp "<<myptcl.species<<"\t"<<plist[i].species<<endl;
        if(plist[i].species == myptcl.species){
            //cout<<"species is "<<i<<endl;
            //cout<<"----------------------------"<<endl;
            plist[i].push(myptcl);
            //cout<<"plist i size now "<<plist[i].member.size()<<endl;
        }
    }
}

void __Cellinfo::pop(const int& sp, const int& offset){
    //cout<<"in pop"<<endl;
    plist[sp].pop(offset);
}



void __Cellinfo::advance(const double& delta_t, vector<vector<__Ptcls> >& templist, const int& myrank)
{

    for(int i = 0; i < plist.size(); i ++)
    {
        //cout<<"prepare to advance ps "<<i<<"\t"<<plist[i].size()<<endl;
        plist[i].advance(delta_t, templist[i], xymin, xymax, myrank);
    }

}

void __Cellinfo::set_field(const __Cellfield& otherfield){
    fields = otherfield;
}






///网格二维数组群，来表征一个box内的cell信息，cell是其内部的member
class __Cellgroup{
public:
    /// 当前group的gridmin
    __Vect2<int> xygridmin;
    /// 当前group的gridmax
    __Vect2<int> xygridmax;
    /// maxid current for ptcls
    int maxid;
    /// 当前group的xymin,左下角坐标
    __Vect2<double> xymin;
    /// 当前group的xymax，右上角坐标
    __Vect2<double> xymax;
    /// 当前group的cell num，在x 和 y 方向
    __Vect2<int> cellnum;
    /// 当前group的cell宽度，dx 和 dy
    __Vect2<double> dxy;
    /// species 是否启用
    vector<Particles> spptcl;
    /// species total
    int species;
    /// species 的 particle per cell
    vector<int> ppc;
    /// 是否关闭em
    int disable_em = 0; 
    /// bd_type
    int bd_type[4];
    /// dt
    double delta_t;
    ///这个位置用来记录在MPI中的位置
    __Vect2<int> position;
    /// cellinfo 数据
    vector<vector<__Cellinfo> > member;
    vector<vector<__Vect3<double> > > jfield;
    vector<vector<double> >rhofield;
    vector<vector<vector<double> > >density_grid;
    vector<vector<vector<double> > >energy_grid;
    __Cellgroup(){}
    __Cellgroup(const __Vect2<int>& xygridmin, const __Vect2<int>& xygridmax, \
            const __Vect2<double>& xymin, const __Vect2<double>& xymax, __Vect2<int>&\
            cellnum, const int& species, const vector<int>& ppc, const vector<Particles>& spptcl, \
            const double& delta_t, const __Vect2<int> position, const int& maxid):xygridmin(xygridmin), \
                                                                xygridmax(xygridmax), \
                                                                xymin(xymin), \
                                                                xymax(xymax), \
                                                                cellnum(cellnum), \
                                                                delta_t(delta_t), \
                                                                position(position), \
                                                                maxid(maxid), species(species), ppc(ppc){
                                                                this->spptcl = spptcl;
                                                                }

    __Cellgroup(const __Cellgroup& other){
        xygridmin = other.xygridmin;
        xygridmax = other.xygridmax;
        xymin = other.xymin;
        xymax = other.xymax;
        cellnum = other.cellnum;
        delta_t = other.delta_t;
        position = other.position;
        species = other.species;
        ppc = other.ppc;
        spptcl = other.spptcl;
        member = other.member;
    }
//    __Cellgroup(){}
    ~__Cellgroup(){}

    /// 初始化，创建cell，没有粒子
    void init();
    /// 压缩当前box的空间
    void squeeze();
    /// 将粒子密度和电流密度插值到网格
    void weight_to_grid(const int& myrank);
    void weight_to_gridb();
    void weight_to_gridc();
    /// 将场数据插值到粒子
    void weight_to_ptcls();
    /// 设置网格位置position处sp类粒子的速度
    void set_velocity(const int& sp, const __Vect2<int>& position, const __Vect3<double>& velocity);
    /// 设置当前box内sp类粒子的速度
    void set_sp_velocity(const int&sp, const __Vect3<double>& velocity);
    /// 初始化当前box内的粒子
    void init_ptcls();
    /// 推动当前cellgroup内粒子运动
    void advance(const double& delta_t, vector<vector<__Ptcls> >& templist, const int& myrank, const int& n);
    /// 输出当前进程rank，当前步n，文件名前缀fname，数据类型sp到文件
    void dump(const int& rank, const int& n, string fname, const int& sp);
    void dump_ptcls(const int& rank, const int& n, string fname, const int& sp);
    /// 使用这个函数来产生事件, 其实只是设置了events_type变量，过程在cell中进行, 0: synchrotron radiation 1: bremsstrahlung radiation, 2: Breit_Wheeler pair, 3: Bethe-Hetiler pair
    friend void monte_carlo(__Cellfield& mygroup, const int& type, const double& delta_t, \
            vector<vector<double> >& tables1d, \
            vector<vector<vector<double> > >& tables2d);
    /// 将rho, j field 清０
    void clear_rhoj();
    /// 合并粒子
    void merge_ptcls(const int& max_ppc, const double& delta_e, const double& delta_v, const int& sp);
    void set_bd(int bd_types[])
    {
        memcpy(bd_type, bd_types, 4 * sizeof(int));
    }
};

void monte_carlo(__Cellgroup& mygroup, const int& type, const double& delta_t, \
        vector<vector<double> >& tables1d, \
        vector<vector<vector<double> > >& tables2d)
{
    //cout<<"size "<<member[50][30].plist[0].member.size()<<endl;
    switch(type)
    {
        case 0:
            {
                for(int i = 1; i < mygroup.member.size() - 1; i ++)
                {
#pragma omp parallel for
                    for(int j = 1; j < mygroup.member[0].size() - 1; j ++)
                    {
                        //cout<<" monte "<<i<<"\t"<<j<<endl;
                        mygroup.member[i][j].synchrotron(delta_t, tables1d, tables2d);
                    }
                }
            }
            break;
        case 1:
            {
                for(int i = 1; i < mygroup.member.size() - 1; i ++)
                {
#pragma omp parallel for
                    for(int j = 1; j < mygroup.member[0].size() - 1; j ++)
                    {
                        mygroup.member[i][j].bremsstrahlung(delta_t, tables1d, tables2d);
                    }
                }
            }
            break;
        case 2:
            {
                for(int i = 1; i < mygroup.member.size() - 1; i ++)
                {
#pragma omp parallel for
                    for(int j = 1; j < mygroup.member[0].size() - 1; j ++)
                    {
                        mygroup.member[i][j].bw_pair(delta_t, tables1d, tables2d);
                    }
                }
            }
            break;
        case 3:
            {
//#pragma omp parallel for
                for(int i = 1; i < mygroup.member.size() - 1; i ++)
                {
#pragma omp parallel for
                    for(int j = 1; j < mygroup.member[0].size() - 1; j ++)
                    {
                        mygroup.member[i][j].bh_pair(delta_t, tables1d, tables2d);
                    }
                }
            }
            break;
        default:
            break;
    }
}


void __Cellgroup::clear_rhoj()
{
    __Vect3<double> zero(0.0, 0.0, 0.0);
    for(int i = 0; i < jfield.size(); i ++)
    {
        for(int j = 0; j < jfield[i].size(); j ++)
        {
            jfield[i][j] = zero;
            rhofield[i][j] = 0.0;
        }
    }
    clear_vec(density_grid);
    for(int i = 0; i < species; i ++)
    {
        density_grid.push_back(rhofield);
    }
}


void __Cellgroup::dump_ptcls(const int& rank, const int& n, string fname, const int& sp)
{
    stringstream tempstream;
    string tempint;
    if(n < 10){
        tempstream<<"0000";
    }
    else if(n < 100){
        tempstream<<"000";
    }else if(n < 1000){
        tempstream<<"00";
    }else if(n < 10000){
        tempstream<<"0";
    }
    tempstream<<n;
    tempstream<<"_rank_";
    tempstream<<rank;
    tempstream>>tempint;
    fname += tempint;
    if(sp >= species)
    {
        cout<<"dump species error: max sp "<<ppc.size()<<endl;
        return;
    }
    string filename;
    if(sp < species)
    {
        filename = fname + "_" + spptcl[sp].name + ".txt";
        ofstream myout;
        myout.open(filename.c_str());
        for(int i = 0; i < member.size(); i ++)
        {
            for(int j = 0; j < member[0].size(); j ++)
            {
                __Ptclslist &mylist = member[i][j].plist[sp];
                for(int k = 0; k < mylist.size(); k ++)
                {
                    __Ptcls& myptcl = mylist.member[k];
                    myout<<myptcl<<endl;
                }
            }
        }    
    }
    
}


 

void __Cellgroup::dump(const int& rank, const int& n, string fname, const int& sp)
{
    stringstream tempstream;
    string tempint;
    if(n < 10){
        tempstream<<"0000";
    }
    else if(n < 100){
        tempstream<<"000";
    }else if(n < 1000){
        tempstream<<"00";
    }else if(n < 10000){
        tempstream<<"0";
    }
    tempstream<<n;
    tempstream<<"_rank_";
    tempstream<<rank;
    tempstream>>tempint;
    fname += tempint;
    int dumpstart = 0, dumpend = 2;
    //cout<<"dumping sp "<<sp<<endl;
    string filename;
    if(sp < species)
    {
        //cout<<"sp "<<spptcl[sp].name<<endl;
        filename = fname + "_density_" + spptcl[sp].name + ".txt";
        ofstream myout;
        myout.open(filename.c_str());
        /***
        for(int i = density_grid[sp][0].size() - dumpend; i > 1; i --){
            for(int j = 2; j <= density_grid[sp].size() - dumpend; j ++){
                myout<<density_grid[sp][j][i]<<"\t";
            }
            myout<<endl;
        }
        ***/
        for(int i = member[0].size() - 1; i > 0; i --){
            for(int j = 1; j <= member.size() - dumpend; j ++){
                myout<<member[j][i].density[sp]<<"\t";
            }
            myout<<endl;
        }
        myout.close();
    }
    else if(sp < species * 2 && sp >= species)
    {
        int thissp = sp - species;
        filename = fname + "_ek_" + spptcl[thissp].name + ".txt";
        ofstream myout;
        myout.open(filename.c_str());
        for(int i = energy_grid[thissp][0].size() - dumpend; i > 1; i --){
            for(int j = 2; j <= energy_grid[thissp].size() - dumpend; j ++){
                myout<<energy_grid[thissp][j][i]<<"\t";
            }
            myout<<endl;
        }
        myout.close();
    }
    else
    {
        cout<<"error sp "<<endl;
    }
/***
    
    switch (sp)
    {
        case 0:
            {
                string filename = fname + "_dens_electron.txt";
                ofstream myout;
                myout.open(filename.c_str());
                for(int i = density_grid[sp][0].size() - dumpend; i > 1; i --){
                    for(int j = 2; j <= density_grid[sp].size() - dumpend; j ++){
                        myout<<density_grid[sp][j][i]<<"\t";
                    }
                    myout<<endl;
                }
                myout.close();
                break;
            }
        case 1:
            {
                string filename = fname + "_dens_ion.txt";
                ofstream myout;
                myout.open(filename.c_str());
                for(int i = density_grid[sp][0].size() - dumpend; i > 1; i --){
                    for(int j = 2; j <= density_grid[sp].size() - dumpend; j ++){
                        myout<<density_grid[sp][j][i]<<"\t";
                    }
                    myout<<endl;
                }
                myout.close();
                break;
            }
        case 2:
            {
                string filename = fname + "_dens_positron.txt";
                ofstream myout;
                myout.open(filename.c_str());
                for(int i = density_grid[sp][0].size() - dumpend; i > 1; i --){
                    for(int j = 2; j <= density_grid[sp].size() - dumpend; j ++){
                        myout<<density_grid[sp][j][i]<<"\t";
                    }
                    myout<<endl;
                }
                myout.close();
                break;
            }
        case 3:
            {
                string filename = fname + "_dens_photon.txt";
                ofstream myout;
                myout.open(filename.c_str());
                for(int i = density_grid[sp][0].size() - dumpend; i > 1; i --){
                    for(int j = 2; j <= density_grid[sp].size() - dumpend; j ++){
                        myout<<density_grid[sp][j][i]<<"\t";
                    }
                    myout<<endl;
                }
                myout.close();
                break;
            }
        case 4:
            {
                int thissp = sp - 4;
                string filename = fname + "_ek_electron.txt";
                ofstream myout;
                myout.open(filename.c_str());
                for(int i = energy_grid[thissp][0].size() - dumpend; i > 1; i --){
                    for(int j = 2; j <= energy_grid[thissp].size() - dumpend; j ++){
                        myout<<energy_grid[thissp][j][i]<<"\t";
                    }
                    myout<<endl;
                }
                myout.close();
                break;
            }
        case 5:
            {
                int thissp = sp - 4;
                string filename = fname + "_ek_ion.txt";
                ofstream myout;
                myout.open(filename.c_str());
                for(int i = energy_grid[thissp][0].size() - dumpend; i > 1; i --){
                    for(int j = 2; j <= energy_grid[thissp].size() - dumpend; j ++){
                        myout<<energy_grid[thissp][j][i]<<"\t";
                    }
                    myout<<endl;
                }
                myout.close();
                break;
            }
        case 6:
            {
                int thissp = sp - 4;
                string filename = fname + "_ek_positron.txt";
                ofstream myout;
                myout.open(filename.c_str());
                for(int i = energy_grid[thissp][0].size() - dumpend; i > 1; i --){
                    for(int j = 2; j <= energy_grid[thissp].size() - dumpend; j ++){
                        myout<<energy_grid[thissp][j][i]<<"\t";
                    }
                    myout<<endl;
                }
                myout.close();
                break;
            }
        default:
            break;
    }

***/
}


void __Cellgroup::squeeze(){
    for(int i = 0; i < member.size(); i++){
        for(int j = 0; j < member[i].size(); j++){
            member[i][j].squeeze();
        }
    }
}


/***
void __Cellgroup::weight_to_gridc(){
    // 先清空cell内场数据, 只需要清空j和rho就可以了
    __Vect2<double> zero2(0, 0);
    __Vect2<double> xhalf(0.5 * dxy.member[0], 0.0);
    __Vect2<double> yhalf(0.0, 0.5 * dxy.member[1]);
    clear_rhoj();
    double area = dxy.member[0] * dxy.member[1];
    double frac = _frac2; //_frac2;
    double idtx = dxy.member[0] / delta_t, idty = dxy.member[1] / delta_t;
    /// 保存临时单个粒子的贡献
    double wtx[5], wty[5], wtxc[5], wtyc[5];
    for(int i = 0; i < member.size(); i ++){
        //cout<<i<<endl;
        for(int j = 0; j < member[i].size(); j ++){
            //cout<<"wt "<<i<<"\t"<<j<<endl;
            __Vect2<double>& cellmin = member[i][j].xymin;
            for(int k = 0; k < member[i][j].plist.size(); k ++){
                //cout<<"sp "<<k<<endl;
                __Ptclslist& mylist = member[i][j].plist[k];
                for(int l = 0; l < mylist.member.size(); l ++)
                {
                    /// 对相邻的三个cell中的场进行插值到粒子
                    /// 这里采用相对坐标local
                    __Ptcls& myptcl = mylist.member[l];
                    /// 两个是坐标是粒子最近格点，另一个是第一次的格点
                    int indx_new = floor((myptcl.position.member[0] + dxy.member[0] * 0.5 - \
                                member[i][j].xymin.member[0]) / dxy.member[0]) + i, \
                    indy_new = floor((myptcl.position.member[1] + dxy.member[1] * 0.5 - \
                                member[i][j].xymin.member[1]) / dxy.member[1]) + j;
                    double cw = myptcl.charge * myptcl.weight / area, wtarea = myptcl.weight / area;
                    __Vect2<double> ratio_new((member[indx_new][indy_new].xymin - myptcl.position) / dxy);
                    __Vect2<double> ratio_newc;
                    ratio_newc.member[0] = (ratio_new.member[0] > 0)?(0.5 - ratio_new.member[0]):\
                                           (0.5 + ratio_new.member[0]);
                    ratio_newc.member[1] = (ratio_new.member[1] > 0)?(0.5 - ratio_new.member[1]):\
                                           (0.5 + ratio_new.member[1]);
                    grid_wxy(ratio_new.member[0], ratio_new.member[1], wtx, wty);
                    grid_wxy(ratio_newc.member[0], ratio_newc.member[1], wtxc, wtyc);
                    
                    double jx = cw * myptcl.velocity.member[0] * frac, jy = cw * myptcl.velocity.member[1] * frac, \
                                 jz = cw * myptcl.velocity.member[2] * frac;
                    //cout<<"wtx1 "<<wtx0[2]<<"\t"<<wtx0[1]<<"\t"<<wtx0[3]<<endl;
                    /// 相邻的3个个点，xy方向都可能有影响
                    int cxmin = _xmin, cxmax = _xmax, cymin = _ymin, cymax = _ymax;
                    for(int s = cxmin; s <= cxmax; s ++)
                    {
                        if(indx_new + s < 1 || indx_new + s > jfield.size() - 1)
                        {
                            continue;
                        }
                        for(int t = cymin; t <= cymax; t++)
                        {
                            if(indy_new + t < 1 || indy_new + t > jfield[0].size() - 1)
                            {
                                continue;
                            }
                            jfield[indx_new + s][indy_new + t].member[0] += jx * wtxc[s + 2] * wty[t + 2];
                            jfield[indx_new + s][indy_new + t].member[1] += jy * wtx[s + 2] * wtyc[t + 2];
                            jfield[indx_new + s][indy_new + t].member[2] += jz * wtxc[s + 2] * wtyc[t + 2];
                            //cout<<" jfield "<<jfield[indx_new + s][indy_new + t]<<endl;

                            rhofield[indx_new + s][indy_new + t] += wtarea * myptcl.charge * wtx[s + 2] * wty[t + 2];
                            density_grid[k][indx_new + s][indy_new + t] += wtarea * wtx[s + 2] * wty[t + 2];
                            energy_grid[k][indx_new + s][indy_new + t] += wtarea * myptcl.energy * \
                                                                          wtxc[s + 2] * wtyc[t + 2];
                        }
                    }
                }
            }
        }
    }
    // 先使用PIC 1st order，这样简单，每个粒子仅仅插值到当前网格的格点上
    // 开始weighting
}

***/
void __Cellgroup::weight_to_grid(const int& myrank){
    // 先清空cell内场数据, 只需要清空j和rho就可以了
    __Vect2<double> zero2(0, 0);
    __Vect2<double> xhalf(0.5 * dxy.member[0], 0.0);
    __Vect2<double> yhalf(0.0, 0.5 * dxy.member[1]);
    clear_rhoj();
    double area = dxy.member[0] * dxy.member[1];
    double frac = _frac2;
    double idtx = 1.0 / dxy.member[0] / delta_t, idty = 1.0 / dxy.member[1] / delta_t, \
                  idxy = 1.0 / dxy.member[0] / dxy.member[1];
    /// 保存临时单个粒子的贡献
    double tempjx[5][5] = {}, tempjy[5][5] = {}, tempjz[5][5] = {};
    double wtx0[5], wty0[5], wtx1[5], wty1[5], wtx0c[5], wty0c[5], wtx1c[5], wty1c[5];

    for(int i = 0; i < cellnum.member[0]; i ++){
        //cout<<i<<endl;
        for(int j = 0; j < cellnum.member[1]; j ++){
            //cout<<"wt "<<i<<"\t"<<j<<endl;
            __Vect2<double>& cellmin = member[i][j].xymin;
            for(int k = 0; k < species; k ++){
                //cout<<"sp "<<k<<endl;
                __Ptclslist& mylist = member[i][j].plist[k];
                if(mylist.size() == 0)
                {
                    continue;
                }
                for(int l = 0; l < mylist.member.size(); l ++)
                {
                    /// 对相邻的三个cell中的场进行插值到粒子
                    /// 这里采用相对坐标local
                    __Ptcls& myptcl = mylist.member[l];
                    /// 两个是坐标是粒子最近格点，另一个是第一次的格点
                    int indx_new = floor((myptcl.position.member[0] + dxy.member[0] * 0.5 - \
                                member[i][j].xymin.member[0]) / dxy.member[0]) + i, \
                    indy_new = floor((myptcl.position.member[1] + dxy.member[1] * 0.5 - \
                                member[i][j].xymin.member[1]) / dxy.member[1]) + j, \
                    indx_old = floor((myptcl.old_pos.member[0] + dxy.member[0] * 0.5 - \
                                member[i][j].xymin.member[0]) / dxy.member[0]) + i, \
                    indy_old = floor((myptcl.old_pos.member[1] + dxy.member[1] * 0.5 - \
                                member[i][j].xymin.member[1]) / dxy.member[1]) + j;
                    indx_old = indx_new, indy_old = indy_new;
                    int dcellx = indx_new - indx_old, dcelly = indy_new - indy_old;
                    double cw = myptcl.charge * myptcl.weight / area, wtarea = myptcl.weight / area;
                    double wx, wy, wz, cwf = myptcl.charge * myptcl.weight * frac;


                    __Vect2<double> localmin_old, localmin_new;
                    localmin_old.member[0] = (indx_old > member.size() - 1)? (cellmin.member[0] + dxy.member[0]):\
                                             (cellmin.member[0]);
                    localmin_old.member[1] = (indy_old > member[0].size() - 1)? (cellmin.member[1] + dxy.member[1]):\
                                             (cellmin.member[1]);
                    localmin_new.member[0] = (indx_new > member.size() - 1)? (cellmin.member[0] + dxy.member[0]):\
                                             (cellmin.member[0]);
                    localmin_new.member[1] = (indy_new > member[0].size() - 1)? (cellmin.member[1] + dxy.member[1]):\
                                             (cellmin.member[1]);
                    __Vect2<double> ratio_old((localmin_old - myptcl.old_pos) / dxy);
                    __Vect2<double> ratio_new((localmin_new - myptcl.position) / dxy);
                    /***
                    __Vect2<double> ratio_oldc;
                    __Vect2<double> ratio_newc;
                    ratio_oldc.member[0] = (ratio_old.member[0] > 0)?(0.5 - ratio_old.member[0]):\
                                           (0.5 + ratio_old.member[0]);
                    ratio_oldc.member[1] = (ratio_old.member[1] > 0)?(0.5 - ratio_old.member[1]):\
                                           (0.5 + ratio_old.member[1]);
                    ratio_newc.member[0] = (ratio_new.member[0] > 0)?(0.5 - ratio_new.member[0]):\
                                           (0.5 + ratio_new.member[0]);
                    ratio_newc.member[1] = (ratio_new.member[1] > 0)?(0.5 - ratio_new.member[1]):\
                                           (0.5 + ratio_new.member[1]);
                    ***/
                    grid_wxy(ratio_old.member[0], ratio_old.member[1], wtx0, wty0);
                    //grid_wxy(ratio_oldc.member[0], ratio_oldc.member[1], wtx0c, wty0c);
                    grid_wxy(ratio_new.member[0], ratio_new.member[1], wtx1, wty1);
                    //grid_wxy(ratio_newc.member[0], ratio_newc.member[1], wtx1c, wty1c);
                    for(int temp = 0; temp < 5; temp ++)
                    {
                        wtx1[temp] -= wtx0[temp];
                        wty1[temp] -= wty0[temp];
                        //wtx1c[temp] -= wtx0c[temp];
                        //wty1c[temp] -= wty0c[temp];
                    }
                    double fjx =  cwf * idty, fjy = cwf * idtx, \
                                 fjz = cwf * idxy * myptcl.velocity.member[2];
                    //cout<<"wtx1 "<<wtx0[2]<<"\t"<<wtx0[1]<<"\t"<<wtx0[3]<<endl;
                    /// 相邻的3个个点，xy方向都可能有影响
                    int cxmin = _xmin + (dcellx - 1) / 2, cxmax = _xmax + (dcellx + 1) / 2, \
                                cymin = _ymin + (dcelly - 1) / 2, cymax = _ymax + (dcellx + 1) / 2;
                    cxmin = _xmin, cymin = _ymin, cxmax = _xmax, cymax = _ymax;
                    for(int s = cxmin; s <= cxmax; s ++)
                    {
                        for(int t = cymin; t <= cymax; t++)
                        {
                            wx = wtx1[s + 2] * (wty0[t + 2] + 0.5 * wty1[t + 2]);
                            wy = wty1[t + 2] * (wtx0[s + 2] + 0.5 * wtx1[s + 2]);
                            wz = wtx0[s + 2] * (wty0[t + 2] + 0.5 * wty1[t + 2]) + \
                                 wtx1[s + 2] * (0.5 * wty0[t + 2] + 1.0 / 3 * wty1[t + 2]);
                            tempjx[s + 2][t + 2] = tempjx[s + 1][t + 2] - fjx * wx;
                            tempjy[s + 2][t + 2] = tempjy[s + 2][t + 1] - fjy * wy;
                            tempjz[s + 2][t + 2] = fjz * wz;
                            //cout<<"jx ,jy, jz "<<tempjx[s + 2][t + 2]<<"\t"<<tempjx[s + 2][t + 2]<<"\t"<<\
                                tempjz[s + 2][t + 2]<<endl;
                        }
                    }
                    for(int s = cxmin; s <= cxmax; s ++)
                    {
                        if(indx_old + s < 0 || indx_old + s > jfield.size() - 1)
                        {
                            continue;
                        }
                        for(int t = cymin; t <= cymax; t++)
                        {
                            if(indy_old + t < 0 || indy_old + t > jfield[0].size() - 1)
                            {
                                continue;
                            }
                            //cout<<"indx , y "<<indx_old + s<<"\t"<<indy_old + t<<endl;
                            jfield[indx_old + s][indy_old + t].member[0] += tempjx[s + 2][t + 2];
                            jfield[indx_old + s][indy_old + t].member[1] += tempjy[s + 2][t + 2];
                            jfield[indx_old + s][indy_old + t].member[2] += tempjz[s + 2][t + 2];
                            //cout<<" jfield "<<jfield[indx_new + s][indy_new + t]<<endl;
                            rhofield[indx_old + s][indy_old + t] += wtarea * myptcl.charge * wtx0[s + 2] * wty0[t + 2];
                            density_grid[k][indx_old + s][indy_old + t] += wtarea * wtx0[s + 2] * wty0[t + 2];
                            energy_grid[k][indx_old + s][indy_old + t] += wtarea * myptcl.energy * \
                                                                            wtx0[s + 2] * wty0[t + 2];
                        }
                    }
                }
            }
        }
    }
       // 先使用PIC 1st order，这样简单，每个粒子仅仅插值到当前网格的格点上
    // 开始weighting
}


void __Cellgroup::weight_to_gridb(){
    // 先清空cell内场数据, 只需要清空j和rho就可以了
    __Vect2<double> zero2(0, 0);
    __Vect2<double> xhalf(0.5 * dxy.member[0], 0.0);
    __Vect2<double> yhalf(0.0, 0.5 * dxy.member[1]);
    clear_rhoj();
    double area = dxy.member[0] * dxy.member[1];
    for(int i = 1; i < member.size() - 1; i ++){
        //cout<<i<<endl;
        for(int j = 1; j < member[i].size() - 1; j ++){
            //cout<<"wt "<<i<<"\t"<<j<<endl;
            __Vect2<double>& localmin = member[i][j].xymin;
            for(int k = 0; k < member[i][j].plist.size(); k ++){
                //cout<<"sp "<<k<<endl;
                __Ptclslist& mylist = member[i][j].plist[k];
                for(int l = 0; l < mylist.member.size(); l ++)
                {
                    //cout<<"ptcls "<<l<<" size "<<mylist.member.size()<<endl;
                    /// 对相邻的三个cell中的场进行插值到粒子
                    /// 这里采用相对坐标local
                    __Ptcls& myptcl = mylist.member[l];
                    double cw = myptcl.charge * myptcl.weight / area;
                    __Vect2<double>& pcoords = myptcl.position;
                    double jx = cw * myptcl.velocity.member[0], \
                                jy = cw * myptcl.velocity.member[1], \
                                jz = cw * myptcl.velocity.member[2];
                    int indx = floor((myptcl.position.member[0] + dxy.member[0] * 0.5 - \
                                localmin.member[0]) / dxy.member[0]) + i;
                    int indy = floor((myptcl.position.member[1] + dxy.member[1] * 0.5 - \
                                localmin.member[1]) / dxy.member[1]) + j;
                    for(int s = _xmin; s <= _xmax; s ++)
                    {
                        for(int t = _ymin; t <= _ymax; t ++)
                        {
                            __Vect2<double> motion = xhalf * 2.0 * s + yhalf * 2.0 * t;
                            jfield[indx + s][indy + t].member[0] += jx * (weightfunc(pcoords, localmin + xhalf + \
                                    motion, dxy) + weightfunc(myptcl.old_pos, localmin + xhalf + motion, dxy)) * 0.5;
                            jfield[indx + s][indy + t].member[1] += jy * (weightfunc(pcoords, localmin + yhalf + \
                                    motion, dxy) + weightfunc(myptcl.old_pos, localmin + yhalf + motion, dxy)) * 0.5;
                            jfield[indx + s][indy + t].member[2] += jz * (weightfunc(pcoords, localmin + yhalf + xhalf + \
                                    motion, dxy) + weightfunc(myptcl.old_pos, localmin + xhalf + \
                                        yhalf + motion, dxy)) * 0.5;
                            density_grid[k][indx + s][indy + t] += myptcl.weight * weightfunc(pcoords, localmin + \
                                    motion, dxy) / area;
                            energy_grid[k][indx + s][indy + t] += myptcl.weight * myptcl.energy *
                                weightfunc(pcoords, localmin + motion, dxy) / area;
                        }
                    }
                    rhofield[i][j] += cw * weightfunc(pcoords, localmin, dxy);
                    rhofield[i+1][j] += cw * weightfunc(pcoords, localmin + xhalf * 2.0, dxy);
                    rhofield[i][j+1] += cw * weightfunc(pcoords, localmin + yhalf * 2.0, dxy);
                    rhofield[i+1][j+1] += cw * weightfunc(pcoords, localmin + (xhalf + yhalf) * 2.0, dxy);
                }
            }
        }
    }

}


void __Cellgroup::advance(const double& delta_t, vector<vector<__Ptcls> >& templist, const int& myrank, const int& n)
{
    // 计算去掉guardcell
    //cout<<"in group advance, cell size "<<member.size()<<endl;
    if(!disable_em)
    {
        for(int i = 1; i < member.size() - 1; i ++){
            for(int j = 1; j < member[i].size() - 1; j++){
            //cout<<"advanceing cell i, j "<<i<<j<<endl;
                for(int k = 0; k < species; k ++)
                {
                    member[i][j].density[k] = (density_grid[k][i][j] + density_grid[k][i+1][j] + \
                            density_grid[k][i][j+1] + density_grid[k][i+1][j+1]) * 0.25;
                }
            }
        }
    }
   for(int i = 1; i < cellnum.member[0] - 1; i ++){
        for(int j = 1; j < cellnum.member[1] - 1; j++){
            //cout<<"cell id "<<i<<"\t"<<j<<endl;
            member[i][j].advance(delta_t, templist, myrank);
        }
    }
    
 
    
    
    /***
    if(bd_type[0] == 2)
    {
        for(int i = 0; i < jfield.size(); i ++)
        {
            jfield[i][0] = jfield[i][jfield[0].size() - 2];
        }
    }
    if(bd_type[1] == 2)
    {
        for(int i = 0; i < jfield.size(); i ++)
        {
            jfield[i][jfield[0].size() - 1] = jfield[i][1];
        }
    }
    ****/
            
    //cout<<"------------------------------------"<<templist[0].size()<<endl;
}

void __Cellgroup::set_sp_velocity(const int& sp, const __Vect3<double>& velocity){
    for(int i = 1; i < member.size() - 1; i ++){
        for(int j = 1; j < member[i].size() - 1; j ++){
            member[i][j].set_velocity(sp, velocity);
        }
    }
}


void __Cellgroup::weight_to_ptcls(){
    for(int i = 1; i < member.size() - 1; i ++){
        for(int j = 1; j < member[i].size() - 1; j ++){
            member[i][j].weight_field_to_ptcls();
        }
    }
}





void __Cellgroup::set_velocity(const int& sp, const __Vect2<int>& position, \
        const __Vect3<double>& velocity){
    member[position.member[0]][position.member[1]].set_velocity(sp, velocity);
}

void __Cellgroup::init(){
    // 非均匀网格初始化，怎么初始化？？？
    // 就先采用均匀的吧
    // 这里包含了guard cell i = 0 & i = gridmax， j = 0 & j = max
    // 计算的时候不要包含这些区域就可以了
    __Vect2<double> zeros(0.0, 0.0);
    __Vect2<int> cposition(0, 0);
    if(dxy == zeros)
    {
        __Vect2<double> dxy = xymax - xymin;
        dxy.member[0] = dxy.member[0] / cellnum.member[0];
        dxy.member[1] = dxy.member[1] / cellnum.member[1];
    }
    for(int i = 0; i < cellnum.member[0]; i ++){
        double xbase = xymin.member[0] + dxy.member[0] * double(i);
        // 同一个i共享同一个x
        vector<__Cellinfo> itemp;
        // itemp存储当前i行的初始化cell
        for(int j = 0; j < cellnum.member[1]; j ++){
            double ybase = xymin.member[1] + dxy.member[1] * double(j);
            // xy base 是当前网格的左下角坐标
            // 下面c开头代表cell
            __Vect2<double> cxymin(xbase, ybase);
            __Vect2<double> cxymax(xbase + dxy.member[0], ybase + dxy.member[1]);
            // cposition 是这个group中的相对位置，MPI时使用更加方便
            cposition.member[0] = i;
            cposition.member[1] = j;
            __Cellinfo temp(cposition, cxymin, cxymax, species);
#ifdef __USRDENS
            calc_density(cxymin, temp.density);
#endif
            // 只是形成了cell没有初始化
            itemp.push_back(temp);
        }
        // i行完毕push到group中
        member.push_back(itemp);
        // 清空temp
        clear_vec(itemp);
        //itemp.clear();
        // 其实这步清空我怀疑是不是多余的？
    }
    // 初始化内部电荷与密度场
    __Vect3<double> v3zero(0.0, 0.0, 0.0);
    for(int i = 0; i < cellnum.member[0] + 1; i ++)
    {
        vector<__Vect3<double>> tempfield;
        vector<double> temprho;
        for(int j = 0; j < cellnum.member[1] + 1; j ++)
        {
            tempfield.push_back(v3zero);
            temprho.push_back(0.0);

        }
        jfield.push_back(tempfield);
        rhofield.push_back(temprho);
        clear_vec(tempfield);
        clear_vec(temprho);
    }
    for(int i = 0; i < species; i ++)
    {
        density_grid.push_back(rhofield);
        energy_grid.push_back(rhofield);
    }

    /***
    // guard cell initialization
    int tppc[4] = {0, 0, 0, 0};
    int tid_now[4] = {0, 0, 0, 0};
    for(int i = 0; i < member.size(); i ++)
    {
        //cout<<"up and down bd"<<endl;
        member[i][0].init(tppc, tid_now);
        // 下边界
        member[i][member[0].size() - 1].init(tppc, tid_now);
        // 上边界
    }
    for(int i = 0; i < member[0].size(); i ++)
    {
        //cout<<"left and right bd "<<endl;
        member[0][i].init(tppc, tid_now);
        // 左边界
        member[member.size()-1][i].init(tppc, tid_now);
        // 右边界
    }
    ***/
}

void __Cellgroup::init_ptcls()
{
    cout<<"init ptcls\t"<<ppc.size()<<endl;
    // cell建立完毕，初始化内部粒子
    for(int i = 0; i < member.size(); i ++){
        // 不能回避掉0和max，guardcell
        // 正常初始化，最后不计算就可以了，不然会出错。。。。
        for(int j = 0; j < member[i].size(); j ++){
            // 采用另一个初始化方法，以粒子种类开始
            //cout<<"init ptcls "<<i<<"\t"<<j<<endl;
            member[i][j].init(ppc, spptcl, disable_em);
        }
    }
}







void __Cellgroup::merge_ptcls(const int& max_ppc, const double& delta_e, const double& delta_v, const int& sp)
{
    for(int i = 0; i < member.size(); i ++)
    {
        for(int j = 0; j < member[i].size(); j ++)
        {
            if(member[i][j].plist[sp].member.size() >= max_ppc)
            {
                member[i][j].plist[sp].merge(delta_e, delta_v, max_ppc);
            }
        }
    }
}
/// 拷贝cellinfo数据到cdata中，MPI传输
void copy_to_cdata(__Cellinfo& cellinfo, __Cdata& cdata)
{
    /***
    cdata.data_rho[0] = cellinfo.fields.rhofieldl.member[0];
    cdata.data_rho[1] = cellinfo.fields.rhofieldl.member[1];
    cdata.data_rho[2] = cellinfo.fields.rhofieldr.member[0];
    cdata.data_rho[3] = cellinfo.fields.rhofieldr.member[1];
    cdata.data_j[0] = cellinfo.fields.jzfield;
    cdata.data_j[1] = cellinfo.fields.jxfield.member[0];
    cdata.data_j[2] = cellinfo.fields.jxfield.member[1];
    cdata.data_j[3] = cellinfo.fields.jyfield.member[0];
    cdata.data_j[4] = cellinfo.fields.jyfield.member[1];
    ***/
    cdata.data_e[0] = cellinfo.fields.ezfield;
    cdata.data_e[1] = cellinfo.fields.exfield.member[0];
    cdata.data_e[2] = cellinfo.fields.exfield.member[1];
    cdata.data_e[3] = cellinfo.fields.eyfield.member[0];
    cdata.data_e[4] = cellinfo.fields.eyfield.member[1];
    cdata.data_b[0] = cellinfo.fields.bzfield;
    cdata.data_b[1] = cellinfo.fields.bxfield.member[0];
    cdata.data_b[2] = cellinfo.fields.bxfield.member[1];
    cdata.data_b[3] = cellinfo.fields.byfield.member[0];
    cdata.data_b[4] = cellinfo.fields.byfield.member[1];
}

/// cline是cdata的一个数组
typedef std::vector<__Cdata> __Cline;

#endif


